As a result of generating a system of equations, we got a matrix written in the form of the Kronecker product of two five-diagonal one-dimensional matrices, obtained from B-splines discretizations.

We now need to solve the resulting system of equations.
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} &\frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{20} & \frac{13}{120} & \frac{1}{120} \\ \end{bmatrix} \otimes \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ \vdots \\ u_{k,l} \\ \vdots \\ u_{N_{x-2},N_y} \\ u_{N_{x-1},N_y} \\ u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) B^x_1(x)*B^y_1(y) dxdy \\ \int BITMAP(x,y) B^x_2(x)*B^y_1(y) dxdy \\ \int BITMAP(x,y) B^x_3(x)*B^y_1(y) dxdy \\\vdots \\ \int BITMAP(x,y) B^x_k(x)*B^y_l(y) dxdy \\\vdots \\\int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_{N_y}(y) dxdy \\\int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_{N_y}(y) dxdy \\\int BITMAP(x,y) B^x_{N_x}(x)*B^y_{N_y}(y) dxdy \\ \end{bmatrix} \)
A solution of a system of equations in which the matrix has a Kronecker product structure is possible in a very short time. What does that mean ‘in a very short time’? The computational cost is expressed in the number of operations such as multiplication or addition of numbers necessary to solve the system of equations. In the case of a system of equations where the matrix has a Kronecker product structure, it is possible to solve the system of equations using an algorithm in which the number of operations \( const*N \) where \( N \) is the number of unknowns (number of bitmap approximation coefficients on the mesh, expressed exactly through \( N=N_x*N_y \), where \( N_x,N_y \) are the mesh sizes, while \( const \) stands for some fixed integer number. The algorithm used in this case is called the alternating-directions algorithm. We consider two steps in the solution process. The first step is to take the first of the sub-matrix and arrange the vector of unknowns and the vector of right-hand sides into many sub-sectors, one vector for each column of the calculation grid elements. It was illustrated in Fig. 1, and in the formula below. \( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} w_{1,1} & w_{1,2} & \cdots & w_{1,N_{y-1}} & w_{1,N_y} \\ w_{2,1} & w_{2,2} & \cdots & w_{2,N_{y-1}} & w_{2,N_y} \\ w_{3,1} & w_{3,2} & \cdots & w_{3,N_{y-1}} & w_{3,N_y} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{N_{x-2},1} & w_{N_{x-2},2} & \cdots & w_{N_{x-2},N_{y-1}} & w_{N_{x-2},N_y} \\ w_{N_{x-1},1} & w_{N_{x-1},2} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x-1},N_y} \\ w_{N_x,1} & w_{N_x,2} & \cdots & w_{N_x,N_{y-1}} & w_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) B^x_1(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_1(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_2(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_2(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_3(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_3(x)*B^y_{N_y}(y)dxdy \\ \vdots \\ \int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_{N_y}(y) dxdy \\ \int BITMAP(x,y) B^x_{N_x}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x}(x)*B^y_{N_y}(y)dxdy \\ \end{bmatrix} \)

We have introduced the auxiliary unknowns here \( w* \) which are used for solutions of the first system of equations. We get original unknowns \( u* \) after solving the second system of equations with v, which the right-hand side will be auxiliary unknowns \( w* \).So, we get a system of equations with a five-diagonal matrix, with many right-hand sides. Each sub-sector of each right-hand side corresponds to one column in the element grid, so it has the specified coordinate \( y \), and the coordinate \( x \) varying from 1 to \( N_x \). The unknowns are similarly ordered \( u* \), in which it the second indexes change with the lines, for example \( w_{1,1}, w_{1,2}, ..., w_{1,N_y} \) while in the columns first indexes change. We solve this system of equations (how to do it we will desribe later), we get solutions \( w* \) and we move on to the second step of the alternating-directions solver algorithm. The second step is then to take an analogous second Kronecker product matrix, taking solutions \( w* \) from the first system of equations, arranging them according to the rows of mesh elements (compare Fig. 1 ), arranged in a similar way unknowns \( u* \) and on solving the obtained system of equations with many right-hand sides. \( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\\cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} & u_{2,1} & \cdots & u_{N_{x-1},1} & u_{N_x,1} \\ u_{1,2} & u_{2,2} & \cdots & u_{N_{x-1},2} & u_{N_x,2} \\ u_{1,3} & u_{2,3} & \cdots & u_{N_{x-1},3} & u_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ u_{1,N_{y-2}} & u_{2,N_{y-2}} & \cdots & u_{N_{x-1},N_{y-2}} & u_{N_{x},N_{y-2}} \\ u_{1,N_{y-1}} & u_{2,N_{y-1}} & \cdots & u_{N_{x-1},N_{y-1}} & u_{N_{x},N_{y-1}} \\ u_{1,N_y} & u_{2,N_y} & \cdots & u_{N_{x-1},N_{y}} & u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} w_{1,1} & w_{2,1} & \cdots & w_{N_{x-1},1} & w_{N_x,1} \\ w_{1,2} & w_{2,2} & \cdots & w_{N_{x-1},2} & w_{N_x,2} \\ w_{1,3} & w_{2,3} & \cdots & w_{N_{x-1},3} & w_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{1,N_{y-2}} & w_{2,N_{y-2}} & \cdots & w_{N_{x-1},N_{y-2}} & w_{N_{x},N_{y-2}} \\ w_{1,N_{y-1}} & w_{2,N_{y-1}} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x},N_{y-1}} \\ w_{1,N_y} & w_{2,N_y} & \cdots & w_{N_{x-1},N_{y}} & w_{N_x,N_y} \\ \end{bmatrix} \)

In this second system of equations, every sub-sector, of every right-hand side corresponds to one row in the element grid, so it has a fixed coordinate \( x \), and the coordinate \( y \) changing from 1 to \( N_y \). The unknowns \( u* \) are similarly ordered. There, the first indexes change with lines, for example \( w_{1,1}, w_{2,1}, ..., w_{N_x,1} \) while in the columns the second indexes change. We will solve each of these two systems of equations with many right-hand sides using the Gaussian elimination algorithm for the band matrix which has a linear computational cost.